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SUMMARY 

A  method  for  the  prediction  of  the  response  of  an  elastic  wing  to 
random  loads  at  flight  conditions  using  rigid  model  wind  tunnel  pressure 
fluctuation  measurements  is  presented.  The  wing  is  divided  into  panels  or 
elements,  and  the  load  is  computed  from  measured  pressure  fluctuations  at 
the  centre  of  each  panel.  A  series  representation  with  terms  of  the  correlated 
noise  type  is  used  to  curve  fit  the  experimentally  determined  pressure  power 
spectra.  Two  methods  are  used  to  calculate  the  random  load  spectrum:  a 
constant  correlation  approximation,  and  an  exponential  spatially  decaying 
cross-power  spectrum  model  for  the  pressure.  The  coupling  between  the 
structural  dynamics  and  aerodynamics  of  a  vibrating  wing  is  taken  into 
consideration  using  the  doublet-lattice  method  for  computing  the  unsteady 
aerodynamic  forces.  The  acceleration  and  displacement  response  spectra 
have  been  computed  for  the  F-4E  aircraft  for  various  Mach  numbers, 
dynamic  pressures  and  flight  altitudes.  The  importance  of  the  unsteady  aero¬ 
dynamic  loads  induced  by  the  vibration  of  the  wing  and  input  load  repre¬ 
sentation  is  illustrated  by  comparing  the  theoretical  predictions  with  results 
from  flight 


RESUME 

II  s’agit  d’une  mlthode  pour  predire  la  rgponse  d’une  aile  eiastique 
a  des  charges  allatoires  en  vol  a  partir  des  mesures  de  fluctuations  de 
pression  sur  une  maquette  rigide  en  souffierie.  L’aile  est  divis£e  en  panneaux 
ou  elements,  et  la  charge  est  calcu!4e  a  partir  des  mesures  de  fluctuations  de 
pression  prises  au  centre  de  chaque  panneau.  Une  representation  aerie  des 
termes  du  type  de  bruit  en  correlation  sert  a  tracer  la  courbe  des  spectres 
de  puissance  pression  determines  experimentalement.  Deux  methodes  sont 
employees  pour  calculer  le  spectre  des  charges  aieatoires:  I’approximation 
par  correlation  des  constantes  et  un  module  exponentiel  du  spectre  de  la 
puissance  decroissante  dans  1’espace  pour  la  pression.  La  relation  entre  la 
dynamique  struct urale  et  l’aerodynamique  d’une  aile  en  vibration  est  prise 
en  compte  au  moyen  de  la  methode  de  reseau  des  doublets  pour  calculer 
les  forces  aerodynamiques  variables.  Les  spectres  de  reponse  en  acceleration 
et  en  deplacement  ont  ete  calculus  pour  le  chasseur  F-4E  a  divers  nombres 
de  Mach,  diverses  pressions  dynamiques  et  altitudes  de  vol.  L’importance 
des  charges  aerodynamiques  variables  induitei  par  la  vibration  de  l’aile  et  la 
representation  des  charges  exercees  est  demontree  par  la  comparaison  des 
predictions  theoriques  avec  les  resultats  des  vols  d’essai. 
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AEROELASTIC  RESPONSE  OF  AN  AIRCRAFT  WING  TO  RANDOM  LOADS 


1.0  INTRODUCTION 

In  References  1  and  2  the  dynamic  response  of  a  wing  to  statistical  non-stationary  random 
loads  is  investigated  using  a  time-segmentation  technique.  The  time  history  of  the  applied  load  is 
segmented  into  a  number  of  time  intervals.  In  each  time  segment  the  non-stationary  load  is  repre¬ 
sented  by  the  product  of  a  deterministic  shaping  function  and  a  statistically  stationary  random  func¬ 
tion.  The  total  response  is  obtained  by  summation  of  the  responses  of  the  wing  from  each  time 
segment.  This  approach  is  useful  when  the  aircraft  engages  in  rapid  manoeuvres.  Because  of  the  com¬ 
plexities  of  the  analysis  only  the  structural  response,  without  consideration  of  the  unsteady  aero¬ 
dynamics  around  the  vibrating  wing,  is  analyzed.  It  is  assumed  that  aerodynamic  effects  can  be 
approximately  accounted  for  by  using  a  damping  value  in  the  dynamic  aerolastic  equations  that  is 
determined  either  experimentally  or  from  estimation. 

In  order  to  make  the  problem  of  the  coupling  between  structural  dynamics  and  aero¬ 
dynamics  of  a  vibrating  wing  more  amenable  to  theoretical  analysis,  the  random  loading  on  the  wing 
is  assumed  to  be  statistically  stationary.  This  assumption  can  be  justified  as  long  as  the  flight  condi¬ 
tions  are  unchanged  within  a  reasonably  long  analysis  time,  when  measurements  are  made  of  either 
the  load  or  the  wing  response.  In  wind  tunnel  testing  statistical  stationarity  requires  long  running 
time  compared  to  the  natural  period  of  the  mode  of  vibration  under  investigation.  The  random  loads 
experienced  by  a  wing  of  an  aircraft  arise  from  turbulence  or  flow  separation  at  buffet  conditions. 

The  use  of  statistical  theory  in  the  investigation  of  buffeting  was  pioneered  by  Liepman  (R  f.  3)  who 
examined  the  problem  of  the  lift  force  exerted  on  a  two-dimensional  thin  airfoil  moving  in  turbulent 
air.  Later,  he  extended  the  method  to  wings  of  finite  span  (Ref.  4).  The  analysis  was  generalized  by 
Ribner  (Ref.  5)  using  a  model  of  turbulence  represented  by  the  superposition  of  plane  sinusoidal 
shear  waves  of  all  arientations  and  wavelengths.  The  correlations  between  flight  and  wind  tunnel  tests 
with  predictions  based  on  statistical  approach  had  been  reported  by  various  investigators  (Refs.  6-10). 
Amongst  the  more  recent  studies  in  this  subject  are  those  of  Mullan  and  Lemley  (Ref.  11),  Hwang 
and  Pi  (Refs.  12, 13),  Jones  (Ref.  14)  and  Butler  and  Spavins  (Ref.  15).  All  the  above  studies  consider 
only  the  structural  response  without  including  the  effects  of  the  unsteady  aerodynamics  around  the 
vibrating  wing.  In  report,  the  method  of  Reference  2  is  extended  to  include  the  coupling  between 
the  structural  vibration  and  aerodynamics  of  a  wing  under  random  loading.  The  panel  method  de¬ 
scribed  in  Reference  2,  which  is  a  modification  of  a  method  proposed  by  Schweiker  and  Davis 
(Ref.  16)  for  the  study  of  the  response  of  shells  to  aerodynamic  noise,  is  most  suitable  for  adapting 
the  doublet-lattice  technique  (Refs.  17-19)  in  computing  the  unsteady  aerodynamics. 

2-0  ANALYSIS 

2.1  Dynamic  Aeroelastic  Equations 

Consider  a  Cartesian  co-ordinate  system  x,  y  and  z  fixed  on  the  wing  as  illustrated  in 
Figure  1.  The  displacement  of  the  wing  can  be  expressed  in  terms  of  a  set  of  normal  co-ordinates 
q«(t)  as: 

z(x,y,t)  =  2  0a(x,y)qa(t)  (1) 

a  “  1 


where  4>a  (x,y )  is  the  mode  shape  function  of  the  ath  mode,  and  I  is  the  number  of  modes  required  to 
adequately  represent  z(x,y,t)  in  the  form  of  a  series.  The  dynamic  aeroelastic  equations  governing 
the  response  of  qa(t)  to  a  load  Ca(t)  is  given  in  generalized  co-ordinate  as: 


Maq<,  +  Cacu  +  K.q.  -  eo(t) 


(2) 
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for  the  ath  mode.  The  dots  denotes  differentiation  with  respect  to  time,  Ma ,  Ca  and  are  the 
generalized  mass,  damping  coefficient  and  stiffness  of  the  ath  mode  respectively.  In  the  above  equa¬ 
tion  the  assumption  of  light  damping  is  made,  thus  cross  damping  terms  do  not  appear. 


Define  the  undamped  natural  frequency  cun  as: 

**0t 


(3) 


and  the  damping  ratio  as: 


Using  the  above  two  expressions  Equation  (2)  can  be  rewritten  as  follows: 


(4) 


+  2faCO„a4,  +  =  —  C„(t) 


(5) 


The  generalized  force  can  be  expressed  in  terms  of  the  fluctuating  pressure  p(x,y,t)  as: 


Ca(t)  =  //  <t>a  (x,y  )p(x,y  ,t)dxdy  (6) 

wing 

planform 


where  the  integration  is  taken  over  the  wing  surface.  Let  La(t>>),  Qa(u>)  and  P(x,y,co)  be  the  Fourier 
transform  of  Ca(t),  qa(t)  and  p(x,y,t)  respectively,  and  introduce  a  frequency  response  function 
Ha(<o)  defined  by: 


Ha(<o)  = 


_ 1 _ 

M“  [W"a  "  w2  +i2?<»WWnJ 


(7) 


Equation  (5)  can  be  written  as  follows: 

Qq(o;)  =  Ha(co)La(oj)  (8) 

In  the  term  8a(t)  there  are  two  types  of  forces  acting  on  the  surface  of  the  wing.  Let  the 
component  of  the  load  due  to  buffeting  or  turbulence  be  represented  by  (t)  and  that  from  the 
unsteady  aerodynamic  forces  due  to  vibration  of  the  wing  be  CA  (t),  then 


«„(t)  =  «°(t)  +  K*(t) 

(9) 

Similarly,  the  pressure  can  be  expressed  as 

P(x,y,t)  =  pD(x,y,t)  +  pA(x,y,t) 

(10) 

and  Equation  (10)  becomes 

Q.(«)  “  Ha(a;)(LD(w)  +  LA(a,)) 

(11) 
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where  l£  (oj)  and  L*  (cj)  are  the  Fourier  transforms  of  8^(t)  and  8*  (t)  respectively. 

2.2  Representation  of  External  Load  on  Wing 

From  Equation  (6),  the  generalized  force  due  to  the  random  loading  can  be  expressed  in 
terms  of  the  pressure  pD  (x,y,t)  as  follows: 


C^(t)  =  //  0a(x,y)PD  (x»y.t)dxdy  (12) 

wing 

planform 

The  power  spectral  density  can  be  written  as: 


S  D(co)  *  ////  0a(x1,y,)flia(x2,y2)  S  D  (x,  ,y,  ,  y2  ,u)  dXj  ,dy,  ,dx2  ,dy2 

P 


(13) 


where  S  D(Xj  ,y(  ,x2  ,y2  ,oj)  is  the  cross  power  spectral  density  of  the  pressure  pD(x,y,t)  between 

p 

points  (Xj  ,yj )  and  (x2  ,y2 )  (Fig.  1). 


In  Figure  2  a  panel  representation  of  the  wing  is  shown.  From  Reference  2,  S  D  (co)  can  be 
written  as:  La 


(oj)  = 


/  /  S  D(x1,y1pt2,y2,a;)dAkdAk 
ak  p 


(14) 


where  Ak  is  the  area  of  the  k-th  panel  and  0k  is  an  average  mode  shape  for  the  ‘a’  mode  within  the 
k-th  panel  defined  by 


*!;  =  r~  f  k 

k  ak 


(15) 


When  the  distance  rks  (Fig.  2)  between  the  centres  of  the  k-th  and  8-th  panels  is  large  (Ref.  2),  an 
approximate  form  for  S  D  (cj)  at  very  low  vibration  frequencies  can  be  obtained  by  assuming  the 

Lo 

pressure  at  the  centre  of  a  panel  to  be  representative  of  the  pressure  field  at  all  points  within  that 
panel.  The  pressure  spectrum  at  any  point  is  given  by  that  at  the  centre  of  the  panel  and  the  pressure 
cross  spectrum  between  two  points  in  the  same  panel  is  taken  to  be  equal  to  the  pressure  spectrum  at 
the  panel  centre.  Equation  (14)  can  then  be  simplified  to  the  following: 

V"’  ■  Ji  V"1  (*•) (16) 

where  Sk_  (o»)  is  the  pressure  spectrum  at  the  centre  of  the  k-th  panel.  A  more  accurate  expression 
p 

for  S  D  (<u)  can  be  obtained  by  writing  S  D(Xj  ,y  j  &2  &2 88 
La  P 

« 

-*|y,-y,l 

8  D(xi«yi’x2*y2»a,>  “  8  D(w)  e  e  cosp(x,  -X2)  (17) 
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0} 

where  a  and  fi  are  spatial  decay  coefficients  in  the  x  and  y  directions  respectively,  p  *  — ,  and  Uc  is 

uc 

the  convection  velocity.  Substituting  Equation  (17)  into  Equation  (16)  and  carrying  out  the  integra¬ 
tion  gives  the  following: 


V"' -  ,f.  V"1  (*-‘)2 11 


(18) 


where 


Ik  =  ~  — - ^~~2 — ~  |  p2  -  a2  +  aXk (a2  +  p2 )  +  e_a xk  [(a2  -  p2  )cos  pXk  -  2 opsin  pXkJ  | 


^  (a2  +  p2  ) 


(19) 


and  Xk  and  Yk  are  the  dimensions  of  the  k-th  panel  (Fig.  2). 

2.3  The  Doublet-Lattice  Method 

The  doublet-lattice  method  is  described  in  References  17-19,  and  in  this  report,  only  a 
brief  outline  of  the  method  and  the  equations  needed  in  the  present  analysis  will  be  given. 

The  wing  surface  is  divided  into  small  elements  or  panels  arranged  in  strips  parallel  to  the 
free  stream  and  all  surface  edges,  fold  lines  and  hinge  lines  lie  on  the  panel  boundaries  (Fig.  3).  In 
each  panel  the  steady  flow  is  represented  by  a  horseshoe  vortex  having  the  bound  vortex  of  the  horse¬ 
shoe  system  lying  along  the  quarter-chord  line.  Superimposed  on  the  bound  vortex  are  acceleration 
potential  doublets  of  uniform  strength.  The  control  point  is  located  at  the  three-quarter  chord  point.  TTie 
normalwash  w’  at  panel  ‘r’  non-dimensionalized  with  respect  to  the  free  stream  velocity  U  is  given 
by  the  following: 


w 


=  ?  DtsACp. 


(20) 


where  Drs  is  the  normalwash  factor,  and  the  summation  is  over  all  the  panels.  In  a  co-ordinate  system 

having  the  midpoint  of  the  sending  element  as  origin  and  the  axes  aligned  with  the  sending  element 
(Fig.  4),  the  normalwash  factor  can  be  written  as  (Ref.  20): 

exp[-iu>(x  -  77tanXs)/u]dTj  (21) 

where 

r,  =  [(y  -  rj)2  +  z2]  *  (22) 

Ax,  is  the  chord  of  the  sending  panel  (Fig.  3),  X  is  the  sweep  angle  of  its  quarter-chord,  and  e  is  its 
semi-width.  T(  and  T2  are  the  direction  cosine  functions  given  as: 


T,  -  cos(7,  -  7t) 


(23) 


T2  -  z  [*008(7,  -  7r)  +  (y  -  r})«n(7,  -  7r)] 


(24) 


where  7  and  yr  are  the  dihedral  angles  at  the  sending  and  receiving  points  on  the  respective  panels. 
Kj  and  k2  are  given  in  Appendix  A.  Equation  (21)  can  be  written  in  the  following  form: 
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D0  +  D,  +  D2 

r*  *ri  ri 


(25) 


where  D„  ,  D.  and  D,  are  the  steady,  planar  and  non-plan ar  normalwash  factors  respectively, 

urs  4rs  ^r$ 

given  in  Appendix  B. 

The  normalwash  boundary  conditions  for  the  lifting  surface  elements  are  determined  by 
assuming  the  surfaces  to  move  normal  to  themselves.  Taking  the  substantial  derivative  of  Equation  (1), 
the  normalwash  at  panel  V  is: 


and 


where 


wf  (t) 

Or 


Wr(t)  =  2  wro(t) 

or 


GOC 


(26) 

(27) 

(28) 


c  is  the  average  chord  and  the  '  denotes  non-dimensional  quantities.  The  Fourier  transform  of  Equa¬ 
tion  (27)  gives: 


(go)  = 

a 


(29) 


2.4  Coupling  between  Structural  Dynamics  and  Unsteady  Aerodynamics  for  a  Vibrating  Wing 


The  external  loading  on  the  wing  due  either  to  buffeting  or  turbulence  causes  vibration  of 
the  wing  which  in  turn  generates  unsteady  aerodynamic  forces  on  the  wing  surface.  Figure  5  shows  a 
block  diagram  illustrating  the  procedure  in  calculating  the  response  of  a  wing  to  random  loading  with 

the  aerodynamics  included  in  a  feedback  loop.  The  aerodynamic  load  LA  (go)  can  be  written  as: 

L^(oj)  =  //  <t>a  (x,y)P(x,y,w)dxdy  (30) 

wing 

planform 

where  P(x,y,to)  =  2  P^x.y.GoJQ^Go)  (31) 


P(x,y,Go)  is  the  Fourier  transform  of  pA  (x,y,t)  in  Equation  (10).  Evaluating  Equation  (30)  over  each 
panel  and  summing  overall  panels,  LA  (go)  becomes: 

l£(w)  =  \  pU2  c2  2  2  0*  ACj  A'  Q'(w)  (32) 

a  &  k  *p 

where  '  denotes  non-dimensional  quantities  and  ACp  is  the  pressure  coefficient. 


The  power  spectrum  S  .  (go)  of  the  aerodynamic  load  can  be  written  as: 
L« 

i  f  (*'.)’  acJ„(“,AcpJ<",so; ("K1 


(33) 


where  Sn'  (cj)  is  the  power  spectral  density  of  Q.(cj)  and  *  denotes  the  complex  conjugate. 
In  non-dimensional  form,  Equation  (8)  can  be  re-written  as: 


q;(w)  -  KaH;(o;)  (l;a(w)  + 


where 


1  pU2c 

Ka  "  «  - T 

2  M“Wna 


H' (co)  -  {1 


X  y  y 

The  primed  quantities  for  the  loads  are  nondimensional  with  respect  to  —  pU  c  .  Using  Equa- 

ml 

tion  (32)  an  expression  for  S„  <  (co)  can  be  obtained  from  the  following: 

SQ-(W){l  -  H^(w)Eaa(w)  -  H*'(w)E*a(«)  +  |h;(oo)|2  |Eoa(a;)|2) 

+  |h;(w)|2  X  =  |H;(a,)|2K2SL,D(aj) 


where 


V">  -  K  f  *5 

To  determine  AC*'  (oj)  combine  Equations  (20),  (29)  and  (31)  and  obtain 


—  +  ik%  -  X  Dk,AC*  (w) 


where  0*  replaced  0*  in  Equation  (27).  Hence,  for  known  Mach  number  and  reduced  frequency, 
Dkt  can  be  obtained  from  the  equations  in  the  previous  section  for  a  given  wing  geometry.  The  mode 
shape  factor  can  be  determined  from  structural  calculations.  With  AC  *  known  SQ »  (oj)  can  be 

pa 

obtained  from  Equation  (37),  and  hence  the  displacement  and  acceleration  spectra. 


3.0  RESULTS  AND  DISCUSSIONS 

The  steps  involved  in  computing  the  response  of  the  wing  are  fairly  straightforward  once 
the  necessary  aeroelastic  data  are  given.  As  an  example,  the  response  of  a  F-4E  wing  is  studied  using 
structural  and  experimental  data  supplied  by  Mullans  and  Lemley  (Ref.  11).  Figure  6  is  taken  from 
Reference  11  and  shows  the  choice  of  panel  dimensions  on  the  wing.  A  total  of  eighteen  panels  are 
used  in  the  computations,  and  the  panel  centres  are  marked  numerically  from  1  to  18. 
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The  first  ten  symmetrical  and  antisymmetric  modes  are  considered  in  this  example.  TTie 
natural  frequencies  cjb  ,  generalized  mass  Mq,  average  mode  shape  0*  and  damping  ratio  ta  tor  three 

dynamic  pressures  Q  -  220, 350  and  470  psf  are  tabulated  in  Mullans’  and  Lemley’s  report  (Ref.  11). 
They  also  obtained  buffet  data  from  wind  tunnel  tests  on  a  10%  scale  rigid  three-dimensional  model. 
For  a  clean  wing  with  zero  leading  edge  and  trailing  edge  deflections,  the  rms  pressure  coefficients 
AC  on  the  wing  are  shown  in  Figure  7.  In  this  particular  example,  AC  varies  practically 

prm  >  Prm« 

linearly  along  both  constant  chord  and  constant  span  directions. 


The  power  spectral  densities  of  the  pressure  fluctuations  at  various  points  on  the  wing  for 
Mach  number  0.7  and  angle  of  incidence  12°  are  also  given  in  Reference  11.  The  transducers  are  not 
located  exactly  at  the  panel  centres  and  there  may  be  more  than  one  transducer  per  panel  or  in  some 
cases,  none  at  all.  Figure  8  indicates  the  transducers  (numbers  correspond  to  these  given  in  Ref.  11) 
whose  measured  outputs  are  used  to  represent  the  fluctuating  pressures  on  the  panels  where  they  are 
located. 


Reference  11  also  gives  the  mode  shape  0*  for  the  symmetric  and  antisymmetric  modes. 

To  obtain  — — ,  a  cubic  spline  algorithm  is  used  in  the  present  study.  For  given  Mach  number  and 
a  x 


reduced  frequency,  the  normalwash  factor  Dfcj  in  Equation  (39)  can  be  obtained  using  the  analysis 
of  Section  2.3.  Upon  substitution  in  Equation  (38),  Equation  (37)  can  be  solved  for  SQ'  (cj)  for 

known  power  spectrum  of  the  external  random  load  S  ,  D  (to).  The  spectra  of  the  pressure  fluctua- 

Lo 

tiorwat  each  of  the  18  panels  are  taken  from  Mullans  and  Lemley  for  q  =  470  psf,  M  *  0.7  and  a  *  12°. 
For  other  flow  conditions,  it  is  assumed  that  the  spectra  remain  the  same.  Following  Reference  2  the 
experimentally  determined  power  spectral  density  curves  from  rigid  model  wind  tunnel  tests  are 
expressed  in  series  form  and  the  coefficients  in  the  series  are  obtained  by  a  curve  fitting  procedure. 
Using  proper  scaling  (Ref.  11)  the  pressure  spectra  at  flight  conditions  can  be  obtained  from  wind 
tunnel  tests. 


The  spectrum  of  the  generalized  force  due  to  the  random  load,  (Eq.  (18))  can  be  calculated 
knowing  the  decay  coefficients  a  and  0  and  also  the  convection  velocity  Uc .  Data  for  a  and  0  from 
flight  tests  are  extremely  scarce.  In  order  to  show  qualitatively  the  effect  of  a  decaying  pressure  field 
and  the  difference  in  the  response  spectrum  as  compared  to  that  for  constant  correlation  assumption 
using  Equation  (16),  the  data  from  Coe  et  al.  (Ref.  22)  are  used.  It  is  assumed  in  this  investigation 
that  0  *  a,  and  a  is  determined  from  measurements  of  surface  pressure  fluctuations  in  a  turbulent 
boundary  layer.  The  decay  coefficient  is  expressed  in  the  form 


where  6  is  the  boundary  layer  thickness,  A  and  B  are  constants  and 


(40) 


(41) 


f  is  the  frequency.  From  Reference  22  A  and  B  are  taken  to  be  0.2  and  .3891  respectively  for 
separated  flow,  and  XQ  -  0.001.  To  calculate  the  boundary  layer  thickness,  the  formula  given  by 
Bies  (Ref.  23)  is  used,  that  is, 
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The  Reynolds  Number  R„  is  calculated  based  on  midspan  chord.  Taking  x  to  be  equal  to  the 

cx 

midspan  chord.  The  value  of  6  thus  calculated  when  substituted  into  Equation  (40)  gives  an  approxi¬ 
mate  average  value  for  a  which  is  assumed  to  be  constant  at  all  the  panels.  It  is  a  function  of  Mach 
number,  frequency  and  flight  altitude  only. 


Figures  9  and  10  show  a  comparison  of  the  response  acceleration  power  spectral  density 
from  flight  tests  (Ref.  11)  with  theoretical  predictions  at  M  =  0.78,  Q  =  220  psf  and  M  =  0.82, 

Q  =  350  psf  respectively.  The  location  of  the  accelerometer  is  at  84%  semispan  and  26%  chord.  The 
constant  correlation  assumption  using  Equation  (16)  for  the  input  load  is  applied  to  all  the  panels. 
Two  sets  of  damping  values  are  used :  a  structural  damping  ratio  of  0.05  for  all  modes,  and  damping 
ratios  determined  experimentally  from  flight  tests  (Ref.  11).  The  ten  symmetric  and  antisymmetric 
modes  have  been  considered  in  the  computations.  It  is  seen  from  these  two  figures  that  at  low  fre¬ 
quencies  below  25  Hz  the  theoretical  predictions  are  fairly  good.  However,  at  higher  frequencies,  the 
predicted  values  are  nearly  tenfold  the  experimental  flight  test  results. 

The  effect  of  altitude  and  Mach  number  on  the  acceleration  and  displacement  spectra  at 
constant  dynamic  pressure  Q  are  shown  in  Figures  11  to  14  for  an  assumed  value  of  0.05  for  the 
structural  damping  ratio.  At  Q  =  220  psf,  the  response  at  low  frequencies  does  not  depend  too 
strongly  on  Mach  number  and  altitude,  while  at  the  higher  frequency  range  of  the  spectra  the  effect 
of  Mach  number  can  be  seen.  It  is  of  interest  to  note  from  these  figures  that  the  response  decreases 
for  increasing  Mach  number  and  altitude  at  low  frequencies  around  10  Hz  while  at  higher  frequencies 
the  reverse  is  true. 

Figures  15  to  18  show  the  response  at  two  Mach  numbers  M  ■  0.51  and  0.78  for  various 
values  of  Q  and  altitude.  They  all  show  that  the  response  is  strongly  dependent  on  Q.  Dividing  the 
response  by  Q2 ,  Figures  19  and  20  show  that  the  curves  can  roughly  be  collapsed  into  a  single  curve 
indicating  a  Q2  dependence.  However,  at  M  =  0.78,  a  similar  behaviour  is  not  observed. 

At  constant  altitude  of  15,000  ft.,  the  variation  of  the  response  with  Mach  number  and  Q 
are  shown  in  Figures  23  and  24.  As  expected,  the  response  increases  with  Mach  number  and  Q 
throughout  the  whole  frequency  spectrum.  The  effect  of  varying  the  structural  damping  is  shown  in 
Figures  25  to  28  for  M  =  0.78,  Q  =  220  psf  and  M  =  0.82,  Q  =  350  respectively. 

To  obtain  better  agreement  between  theoretical  and  experimental  results  than  those  shown 
in  Figures  9  and  10  at  the  higher  frequencies,  Equation  (18)  is  used  to  represent  the  input  load 
spectrum  using  values  of  decay  coefficients  derived  from  Equations  (40)  to  (42)  and  a  value  of  the 
convection  velocity  Uc  =  0.80V.  These  results  are  shown  in  Figures  29  to  32.  From  the  acceleration 
spectra  at  M  *  0.78  and  0.82  (Figs.  29  and  31),  it  can  be  seen  that  the  high  frequency  range  of  the 
spectra  has  been  much  improved  and  agrees  well  with  flight  test  data.  However,  at  low  frequencies, 
the  constant  correlation  assumption  seems  to  give  better  agreement.  This  is  probably  due  to  the 
inaccuracies  of  the  decay  coefficients  obtained  by  extrapolation  of  data  obtained  at  much  higher 
frequencies. 


Figures  33  and  34  also  show  comparisons  between  flight  test  data  and  theoretical  predic¬ 
tions.  These  figures  are  prepared  by  combining  Figures  9  and  10  with  Figures  29  and  31  for  M  ”  0.78, 
Q  -  270  psf  and  M  ■  0.82,  Q  ■  350  psf  respectively.  The  predictions  from  Reference  11  which  con¬ 
sidered  only  the  structural  response  are  also  included  to  show  the  effect  of  unsteady  aerodynamics  in 
the  response  calculations.  At  low  frequencies  below  15  Hx,  improvements  over  predictions  which 
considered  only  the  structural  dynamics  of  the  wing  are  quite  noticeable.  However,  at  the  higher 
frequency  range  of  the  spectra,  unsteady  aerodynamic  effects  become  less  important. 


4.0  CONCLUSIONS 


A  panel  method  for  predicting  the  acceleration  and  displacement  response  of  a  wing  to 
random  loading  has  been  developed  taking  into  consideration  the  coupling  between  the  structural 
dynamics  and  aerodynamics  of  a  vibrating  wing.  The  doublet-lattice  method  is  found  to  be  most 
suitable  in  the  present  panelling  scheme  to  compute  the  unsteady  aerodynamic  forces. 

Using  available  structural  data  for  the  F-4E  wing  together  with  rigid  model  wind  tunnel 
fluctuating  pressure  measurements,  the  displacement  and  acceleration  response  have  been  computed 
for  various  Mach  numbers,  dynamic  pressures  and  flight  altitudes.  The  results  indicate  that  the 
response  is  most  affected  by  changes  in  dynamic  pressure  and  Mach  number. 

At  low  frequencies,  improvements  in  the  response  predictions  over  those  which  considered 
only  the  structural  dynamics  but  excluded  the  unsteady  aerodynamic  forces  generated  by  the 
vibrating  wing  are  quite  significant.  However,  at  the  higher  frequency  range  of  the  spectra,  unsteady 
aerodynamic  effects  become  less  important. 

Two  methods  have  been  used  to  represent  the  spectrum  of  the  input  random  load.  Using 
the  constant  correlation  assumption,  comparison  of  the  computed  acceleration  power  spectra  with 
those  from  flight  tests  shows  that  at  low  frequencies  the  agreement  is  fairly  good,  but  at  high  fre¬ 
quencies,  the  predicted  values  are  nearly  tenfolds  the  experimental  results.  Much  better  agreement 
at  the  higher  frequencies  is  obtained  when  the  exponential  spatially  decaying  form  of  the  cross-power 
spectral  density  for  the  fluctuating  pressures  is  used.  However,  this  results  in  discrepancies  at  low 
frequencies  which  are  probably  due  to  the  extrapolation  of  the  decay  coefficients  from  data  obtained 
at  much  higher  frequencies.  The  results  illustrate  the  importance  of  a  proper  representation  of  the 
input  load  spectrum. 
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FIG.  3:  REPRESENTATION  OF  WING  BY  PANELS  AND  LOCATION  OF  VORTICES. 
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FIG.  4:  WING  AND  ELEMENT  CO-ORDINATES 
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FIG.  5:  BLOCK  DIAGRAM  FOR  WING  RESPONSE  TO  RANDOM  LOADING 
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FIG.  9:  COMPARISON  OF  RESPONSE  ACCELERATION  POWER  SPECTRAL  DENSITY 
FROM  FLIGHT  TESTS  WITH  THEORETICAL  PREDICTIONS  USING  CONSTANT 
CORRELATION  ASSUMPTION  AT  M  «  0.78  AND  Q  -  220  PSF 
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FIG.  11:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  WITH 
MACH  NUMBER  AND  FLIGHT  ALTITUDE  FOR  Q  -  220  PSF 
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FIG.  12:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
MACH  NUMBER  AND  FLIGHT  ALTITUDE  FOR  Q  »  220  PSF 
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FIG.  13:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  WITH 
MACH  NUMBER  AND  FLIGHT  ALTITUDE  FOR  Q  -  350  PSF 


DISPLACEMENT  POWER  SPECTRAL  DENSITY,  FT*»2/H2 


-26- 


FREQUENCY,  HZ 


FIG.  14:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
MACH  NUMBER  AND  FLIGHT  ALTITUDE  FOR  Q  -  350  PSF 
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FIG.  16:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  WITH 
DYNAMIC  PRESSURE  AND  FLIGHT  ALTITUDE  FOR  M  -  0.51 
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FIG.  16:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
DYNAMIC  PRESSURE  AND  FLIGHT  ALTITUDE  FOR  M  -  0.51 
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FIG.  17:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  WITH 
DYNAMIC  PRESSURE  AND  FLIGHT  ALTITUDE  FOR  M  -  0.78 
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FIG.  18:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
DYNAMIC  PRESSURE  AND  FLIGHT  ALTITUDE  FOR  M  =  0.78 
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FIG.  19:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  NORMALIZED 
WITH  RESPECT  TO  Q2  WITH  DYNAMIC  PRESSURE  AND 
FLIGHT  ALTITUDE  FOR  M  -  0.51 
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FIG.  20:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  NORMALIZED 
WITH  RESPECT  TO  Q2  WITH  DYNAMIC  PRESSURE  AND 
FLIGHT  ALTITUDE  FOR  M  =  0.51 
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FIG.  21:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  NORMALIZED 
WITH  RESPECT  TO  Q*  WITH  DYNAMIC  PRESSURE  AND 
FLIGHT  ALTITUDE  FOR  M  -  0.78 
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FIG.  22:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  NORMALIZED 
WITH  RESPECT  TO  Q*  WITH  DYNAMIC  PRESSURE  AND 
FLIGHT  ALTITUDE  FOR  M  -  0.78 
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FIG.  24:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
MACH  NUMBER  AND  DYNAMIC  PRESSURE  AT  15,000  FT.  ALTITUDE 
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FIG.  26:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
STRUCTURAL  DAMPING  FOR  M  -  0.78,  Q  -  220  PSF 
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FIG.  27:  VARIATION  OF  ACCELERATION  POWER  SPECTRAL  DENSITY  WITH 
STRUCTURAL  DAMPING  FOR  M  -  0.82,  Q  ■  350  PSF 
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FIG.  28:  VARIATION  OF  DISPLACEMENT  POWER  SPECTRAL  DENSITY  WITH 
STRUCTURAL  DAMPING  FOR  M  -  0.82,  Q  -  350  PSF 
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FIG.  29:  COMPARISON  OF  RESPONSE  ACCELERATION  POWER  SPECTRAL  OENSITY 
FROM  FLIGHT  TESTS  WITH  THEORETICAL  PREDICTIONS  USING  EXPONENTIAL 
SPATIALLY  DECAYING  PRESSURE  CROSS-SPECTRUM  ASSUMPTION  FOR 

M  -  0./8,  Q  *  220  PF 
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FIG.  30:  DISPLACEMENT  POWER  SPECTRAL  DENSITY  USING  EXPONENTIAL 
SPATIALLY  DECAYING  PRESSURE  CROSS-SPECTRUM  ASSUMPTION  FOR 

M  -  0.78,  Q  -  220  PF 
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FIG.  31:  COMPARISON  OF  RESPONSE  ACCELERATION  POWER  SPECTRAL  DENSITY 
FROM  FLIGHT  TESTS  WITH  THEORETICAL  PREDICTIONS  USING  EXPONENTIAL 
SPATIALLY  DECAYING  PRESSURE  CROSS-SPECTRUM  ASSUMPTION  FOR 

M  -  0.82,  Q  *  350  PSF 
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FIG.  32:  DISPLACEMENT  POWER  SPECTRAL  DENSITY  USING  EXPONENTIAL 
SPATIALLY  DECAYING  PRESSURE  CROSS-SPECTRUM  ASSUMPTION  FOR 

M  -  0.82.  Q  -  350  PSF 
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FIG.  33:  COMPARISON  OF  RESPONSE  ACCELERATION  POWER  SPECTRAL  DENSITY 
FROM  FLIGHT  TESTS  WITH  THEORETICAL  PREDICTIONS  FOR 
M  *  0.78,  Q  *  220  PSF 


FIG.  34:  COMPARISON  OF  RESPONSE  ACCELERATION  POWER  SPECTRAL  DENSITY 
FROM  FLIGHT  TESTS  WITH  THEORETICAL  PREDICTIONS  FOR 
M  -  0.82,  Q  -  350  PSF 
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APPENDIX  A 


EXPRESSIONS  FOR  K,  AND  K2  IN  THE  EQUATION  FOR  NORMALWASH 


In  Equation  (21),  the  terms  K  and  K,  are  given  by  Landahl  (Ref.  20)  as  follows  (with 
sign  reversed): 


K,  "  'I, 


exp(-ikj  u,  )Mr, 


0“ 


K,  =  31,  + 


ikj  M2r^  exp(-ik,  Uj )  Mr! 


R2  (l  +u2)* 


-r  M) 


MrjU,  expt-ikjUj) 


where 


0  -  (1-  M2) 


R  =  [(x  -  fjtan\s)2  +  02r2J  ' 


U,  =  (MR  -  x  +  T}tanXs)/02r, 


k.  =  - 

1  U 


For  u,  >  0  (Ref.  19), 


MU1  )  =  1 


("4 


-  ikjlpfuj.k,)  expf-ikjU,) 


where 


u  a  (nc-ik,)exp(-ncu, ) 
V“i  - - — : - - 


n2c2  +k2 
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0.2418619 


-2.7918027 


24.991079 


=  -111.59196 


271.43549 


=  -305.75288 


=  -41.1863630 


=  545.98537 


=  -644.78155 


aio 


328.72755 


64.279511 


0.372 


*i’ki)  “  |(2  +  ik,u, )  A  ^  ikj  I0(Uj . 

I  \  (i-i) 


+  k, J0(u,  ,k,  )>exp(-ik,u, ) 


(A10) 


where 


J0(ui  »ki ) 


,,  ,,  , 

"'(„v*kj)!  '  v  ' 


ik,  ^2nc  +  Uj  ^n2c2  +ki)j| 


For  u,  <  0, 


Ij  (Uj  ,kj )  =  2118  1,(0*,)  -  R8I,  (~u,  ,k, )  +  ilm  I,  (-u,,k,) 


(All) 


(A12) 


I2(u,  Jk, )  =  2R8  l2(0,k, )  -  R8I2  (-u,,k,)  +  ilm I,  (-u,,k,) 


(A13) 


For  steady  states,  K,  and  K2  are  given  as: 


K,n  =  -1  -  (x  -rjtanX  )/R 


K20  *  2  +  (x  -rjtanXj)  12  + 


P2rl\ 

V } 


(A14) 


(A15) 


APPENDIX  B 


EXPRESSIONS  FOR  NORMALWASH  FACTOR 


In  Equation  (25),  the  normalwash  factor  is  given  as: 


Drs  =  D0  +  D1  +  °2 
rs  rs  *rs 


From  Reference  (21),  the  steady  normalwash  factor  Dn  can  be  derived  as  (Fig.  Bl): 

urs 

D0  =  —  (Vysin7r  -  Vzcos7r) 


where 


Vv  =  (dbxcosAsin7s  -  dh7sinA) 


cos6b  -  costpb  1  -  cos*pi 


+  1 


s  bz 


R:„  ~ 


4*<  “  4"(<  tR0z) 


and 

cos0b  -  cosy3b 

V  =  (db  sinA- dbxcosAcos7j)  - - - 


1  -  COS*p.  COS0Q  +  1 

4-«  *r?,)  Ri>  4*«+Ry 


The  components  of  db ,  R.  and  R0  are : 

d.  =  Riv  -  R.cosfl.  sinA 

DX  IX  ID 

dby  =  Riy  "  Ric0S®bC08Ac0S7s 

dbz  =  Rjz  -  RicosflbcosAsin7s 

and 

Rjx  =  xo  +  etanA 

Riy  =  Fo  +  ec°8T$ 

Riz  "  zo  +  e8in')r1 
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where 


X„  =  — 

0  P 


tan  A  =  l/0tanX 


P  =  yll-  M2 


Finally,  the  angle  relations  are  given  as: 


(B13) 


(B14) 


(B15) 


cosd.  =  -  (R.  sinA  +  R.  cosAcos7  +  R.  cosAsury  ) 


R;  1X 


(B16) 


cos<pb  =  —  (R0xsinA  +  RQycosAcos7s  +  RQzcosAsin7s) 


(B17) 


COS^j  =  - 


(B18) 


.  Ox 

cos6n  ~  - 

0  Ro 


From  Reference  (19),  the  planar  normalwash  factor  D,  is  given  as: 

1  rs 

The  terms  A, ,  Bj  and  Cj  are  determined  from  the  following: 

Al  =  -V{P,(-«)-2P1(0)  +  PI(e)} 

2e2 

B,  -  ^  {P,(e)- P,(-e)) 


(B19) 


+  2eA. >  (B20) 


(B21) 


(B22) 


C,  =  P,(0) 


(B23) 


where  P,(e),  Pt(-e)and  P(0)  are  the  outboard,  inboard  and  center  values  of  P(rj)  defined  as: 


’,(»?)  -  T,  Jk,  exp  ^  (x-TjtanXs)J  -  K,0j 


(B24) 


which  is  obtained  by  subtracting  the  steady  value  K|0  from  the  numerator  of  the  first  term  of  Equa¬ 
tion  (21).  The  term  F  in  Equation  (B20)  is  given  by: 


where  the  principal  value  of  the  arctangent  is  taken  in  the  range  (0,ir).  When  z  =  0, 


and  for  z  «  l,y2  +  z2  >  e2, 


and 


(B26 


(B27 


(B28 


The  non-planar  normal  wash  factor  D2  is: 

^rs 


+  (y2  +  z2  -  ye)B2 


+  (y-  e)C 


(B29 


A2 ,  B2  and  C,  are  determined  from  Equations  (B21)  to  (B24)  with  the  subscripts  in  those  equations 
changed  to  *2  .  For  small  values  of  z,  the  following  equation  for  D,  is  used: 


Ax,e 

2(y2 

+  z2  +  e2 )  (e2  A 

2  +C2)  +  4ye2B2 

8x(y2  +z2  -  e2) 

[(y  +  e)2  +  Z2] 

[(y  -  e)2  +  z2] 

+  z2  )A2  +  yB2  + 


(B30) 


where  a  is  given  by  Equation  (B28).  This  equation  is  generally  used  except  when  (y2  +  z2  -  e2)/ 
2ez  <  0.1  in  which  case  Equation  (B29)  is  used  instead. 
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FIG.  B1:  HORSESHOE  VORTEX  SYSTEM 
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